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1 Introduction 


Open region scattering problems ave frequently modelled using an integral 
equation formulation[l] and solved using the method of moments. For an 
inhomogeneous dielectric scatterer, the integral equation must be discretized 
over the entire volume of the object. Since this leads to a full matrix system, 
the integral equation technique proves to be expensive in terms of storage 
and solution time. Differential equation methods, such as finite elements, 
have therefore been tried as an alternative solution approach. Although 
these techniques give rise to matrix systems that are considerably larger than 
those generated using the integral equation formulation, the finite element 
matrices are highly sparse leading to smaller storage and lower solution time. 
These methods, however, can only solve bounded field problems whereas 
most electromagnetic problems are open boundary-infinite domain types. To 
solve for unbounded field problems, the finite element mesh needs to be 
truncated artificially at some distance from the scatterer. A wide variety of 
methods like ballooning, harmonic series expansions and infinitesimal scaling 
have been developed to accomplish this. However, all of the above methods 
are ‘global’ techniques which produce fully populated submatrices that are 
computationally expensive and spoil the sparse, banded structure of the finite 
element system. 

The most promising method of mesh termination developed so far has 
been the application of an absorbing boundary condition (ABC) on an artifi- 
cial outer boundary to minimise the non-physical reflections from the bound- 
ary. The essence of the method is to force the field components at the mesh 
termination plane to satisfy the differential equation (absorbing boundary 
condition) for an outward travelling wave. The degree of accuracy of the fi- 
nite approximation to the differential equation determines the quality of the 
mesh termination. Stability considerations limit the degree of approxima- 
tion. The advantage of the ABCs over global methods is that they are local 
boundary conditions and hence retain the sparse structure of the finite ele- 
ment formulation. Moreover, the additional computational effort when using 
ABCs is small when compared to a bounded field problem. One disadvan- 
tage, however, is that ABCs are approximate and do not model the exterior 
field exactly. The local boundary conditions should also give rise to inac- 
curacies when travelling waves are excited on the scatterer and substantial 
global coupling exists between widely separated parts of the scatterer. 

In this report, we consider a FE-ABC solution of the scattering by arbi- 
trary three-dimensional structures. The computational domain is discretized 
using edge-based tetrahedral elements. In contrast to the node-based ele- 
ments, edge elements can treat geometries with sharp edges, are divergence- 
less and easily satisfy the field continuity condition across dielectric interfaces 
[1]. They do, however, lead to a higher unknown count but this is balanced 
by the greater sparsity of the resulting finite element matrix. Thus the com- 
putation time required to solve such a system iteratively with a given degree 
of accuracy is less than the traditional node-based approach [2]. 
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The purpose of this report is to examine the derivation and performance 
of the ABCs when applied to two and three dimensional problems and to dis- 
cuss the specifics of our FE-ABC implementation. The two dimensional ones 
investigated here are the well-known Bayliss-Turkel[3] and Engquist-Majda[4] 
absorbing boundary conditions. The three dimensional ABCs presented here 
are those derived by Mittra[5] and Kanellopoulos and Webb [6]. All the ABCs 
discussed here are derived from the Wilcox expansion and can only be ap- 
plied on a spherical outer boundary. Some results are then presented which 
demonstrate that remarkably accurate solutions can be obtained by enforc- 
ing the ABC a small fraction of a wavelength from the scatterer. This is in 
contrast to our experience with two-dimensional applications and is probably 
due to the faster (1 fr rather than 1 !\fr) decay of the propagating fields. 


2 Derivation of the ABCs 

2.1 Two-dimensional ABCs 

Absorbing boundary conditions for two-dimensional problems have been ex- 
tensively studied over the years. The most notable ABC to come out of this 
research has been the one proposed by Bayliss, Gunzburger and Turkel[3] 
who employed an asymptotic analysis to derive a series of local operators. 
Using the pseudo-differential operator theory, Engquist and Majda[4] have 
generated a set of different operators for minimisation of the reflections from 
the mesh termination boundary. However, since the Engquist-Majda ABCs 
are not as accurate as the Bayliss-Gunzburger-Turkel(BGT) ABCs, only the 
BGT operators will be discussed here. 

Consider a perfectly conducting cylindrical scatterer, shown in Figure 1, 
whose cross-section is defined by the contour Ti. Let the exterior region of 
the scatterer lie in the domain Cl. For a TM-polarized incident wave, we need 
to solve the wave equation 


V 2 u + k 2 u = 0 (1) 

in Cl. The wave function u is proportional to the z-component of the scattered 
electric field and satisfies the boundary condition 

u' + ti = 0 onTi (2) 

where u' is the incident wave function. 

Following Wilcox, an asymptotic expression for u can be written as fol- 
lows: 


u(p,<f>) 


,-jkp 


Vp L 




e~ lkp A QnW 
P 1/2 ho P n 



( 3 ) 
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Defining u p as the derivative of u with respect to p, we have from (3) 


Up + jku 


e -jkp 

v 71 


ac{<t>) + 3 


a i{4>) 

p 




( 4 ) 


Therefore, from (4), we obtain 


Up + j ku 



Thus if we neglect terms of the order 0(p ~ 3 / 2 ) and smaller, we obtain u p + 
jku = 0, which is equivalent to the Sommerfeld radiation condition for the 
wave function u in two dimensions. It is shown next that a higher-order 
boundary operator B\ can be obtained that yield terms of the order 0(p ~ 5 / 2 ) 
when applied to u. The operator B\ is given by 

B ‘ = T P +jk + Tp (5) 

It is readily found that, for a given p, B\ introduces a higher order error 
in p -1 than does the Sommerfeld radiation condition such that terms of the 
order 0(p ~ 5//2 ) and smaller are neglected. Continuing along similar lines, the 
next higher order operator B 2 can be derived by first defining v = B\U, and 
then showing that 

{r P +jk+ ^) v = 0{P ' W) 

The second-order absorbing boundary operator B 2 is thus given by 


B 2 “ {r P +jk+ T P , 


9 -i 1 

o' + ]k+ — 

,op 2 p / 


(6) 


In fact, it is shown in [3] that a generalised operator B m can be constructed 
by repeating the above procedure such that 


B„ 


= n 

p=i 


d 

T P +lk + 


2 p - 3/2' 


( 7 ) 


and 


B m u 


O 


1 

p 2 m +\/2 


Since the error between the exact field and the approximated value decreases 
as we employ higher order boundary conditions, the mesh termination bound- 
ary can be brought closer to the scatterer resulting in a smaller number of 
unknowns. However, the higher order conditions make the system of equa- 
tions progressively more ill-conditioned and the choice of the order of the 
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ABC is often limited by this constraint. The second order ABC is usually 
employed in practice. 

An alternate derivation[5] of the boundary operators can also be carried 
out by imposing the requirement that u(p, <f > ) in (3) satisfy the wave equation 
and then deriving a recursion relationship in terms of the angular derivatives. 
The second order ABC thus reduces to 


u p 


( jk 2p + Sjkp* + 8Pp 3 ) U + 
a(p)u + P(p)um 




Ufifi + O 



Higher order boundary operators in the above formulation involve higher 
order angular derivatives and are thus simple to implement numerically. 


2.2 Three dimensional ABCs 

Three dimensional absorbing boundary conditions can be classified into two 
categories: scalar ABCs and vector ABCs. Most of these boundary con- 
ditions are untested till date; it is, therefore, difficult to comment on the 
performance of these boundary operators. However, in tests carried out us- 
ing the vector ABC derived in [6], reliable results have been obtained on 
placing the mesh truncation boundary a fraction of a wavelength distant 
from the scatterer. 

The first section presents a scalar ABC proposed in [5] and the second 
section is devoted to the formulation of vector ABCs derived in [5] and [6]. 


2.2.1 Scalar ABC 


Let u(r,0,4>) satisfy the scalar wave equation (1) and be expressed asymp- 
totically as 


u = 


e~ jkr a n (0, <t>) 

nr* Tl 


(9) 


n=0 


It can then be shown that the coefficients in the series expansion, a n , satisfy 
the following recursion relationship 


—2 jkna n — n(n — l)G n _! + Da n - x 
where D is Beltrami’s operator and is given by 


D = 


1 d 
sin 6 dO 


sin 0 


d 

dO 


+ 


1 d 2 

sin 2 6 d(f> 2 


( 10 ) 


( 11 ) 


Differentiating (9) with respect to r and incorporating the the recursion 
relationship repeatedly in the resulting expression yields 

U r = —jk{a(r)u + /3(r)Du) (12) 

up to and including terms of order r~ 4 . In (12), 


a(r ) = 1 + - 


1 


jkr 


P(r) = 


1 


2 (kr) 2 a(r) 
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2.2.2 Vector ABC 


As for the scalar case, we begin with the Wilcox representation of the electric 
field which has an expansion of 


E(r) = 


r n=o r 


where r,9,<f> are spherical coordinates. From (13), we get 

f / 1\ e- jlr “ nA„, 

V*E={-(^ + -)rx + 7 i}E-— 


(13) 


(14) 


where A nt = r x A n is the transverse component of A„ and, for a vector F, 
DiF is given by 


Di F = 


1 

sinO 

1 

sinO 


d . . 8F*1 


\ dfT ■ 

— — sindr v 

0 + 

i 

; 

i 

<c> 

1 

89 


L w\ 




(15) 


Using the recursion relation 


2jknA n t — n.(n l)A n _i,t -I- D 4 A n —\ 


where 


D 4 A n 


Dg A n 

D</> A n 


(B a; + D, A n ) « + (B <( + A„) <£ 

89 sin 2 9 n sin 2 9 89 

J2_8K 1_ + 2cos9 8A e n 

sin9 84 > sin 2 9 n sm 2 # 84> 


and D is Beltrami’s operator defined in (11), we can derive the representation 
correct to r~ 4 . Applying the recursion relation in ( 14) yields the desired 
relationship for the vector ABC: 


V x E = a(r) E + /3(r)D 4 E 


(16) 


where 


a(r) 

P( r ) 




1 


1 


2 jkr 2 (1 — 1/jkr ) 


An alternate derivation of the three-dimensional vector ABC is given in [6] 
which has a higher order of accuracy and is easier to implement than the 
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previous one. The basic building block is the differential operator Ln, defined 
as 

Ljv( u) = rxVxu-|jH-j u, TV = 0,1,2,... (17) 

where r is a unit vector in the radial direction. It can be shown that for 
TV > 0 and n > 0, 


Ln 



(" - ( 18 > 

+ (19) 


where g is e -jfcr and the subscripts t and r denote the transverse and radial 
components of a vector, respectively. In both cases, L n has the effect of 
multiplying by 1/r while leaving the transverse dependence unchanged. 

The operators Bn, TV = 1,2,... can now be defined such that 

S »(E) = o(^fr) 

where 


Bn g 


A „(M)' 


r-n+1 


(" + 1 - N )( n + 2 - N) ■ ■ ■ (n)g A *^p - 

+ S (n + 1 - N)(n + 2 - N) ■ ■ ■ (n - 1)V, 

( 20 ) 


Since A 0r is zero, it can be seen that the RHS of (20) vanishes for n = 
0,1,..., TV — 1 , i.e., Bn annihilates the first TV terms of the vector expansion 
(13). Thus, Bn = 0 is an approximate absorbing boundary condition on 
the surface 5" of a sphere of radius r. After a certain amount of algebra 
and making use of the fact that E satisfies the vector wave equation, the 
operators B\ and B% can be written as 


Z?i(E) = r x V x E - qE ( + (s - l)V t E r (21) 

B 2 { E) = J(rxVxE) + ^E f + Vx[r(VxE r ] 

+(s - l)V t (V • E t ) + (2 - s)aV t E r (22) 

where a = jk and f5 = 1/(2 jk + 2/r). The operator B 2 can be incorporated 
into a weighted residual formulation for s = 1 but leads to an unsymmetric 
matrix problem. However, for s = 2, the functional yields a sparse, sym- 
metric matrix. Thus, the second order boundary operator can be expressed 
as 


B 2 { E) = -(? x V x E) + P 2 (E) 
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where 


P 2 { E) = aE t + /3Vx [f(V x E) r ] + (s - 1)/?V,(V • E t ) + (2 - s)a(3V t E r ( 23) 

which can be easily incorporated into the surface integral of the electric field 
functional. 


3 Formulation 


3.1 Two-dimensional case 


The scattered field must satisfy the wave equation (1) within the domain 
of interest. The first step in the finite element formulation is to multiply 
(1) with a test function v and integrate the product over the entire problem 
domain. This gives 


J (uV 2 u + k 2 vu ) dS = 0 (24) 

where u consists of a set of known basis functions weighted by its correspond- 
ing unknown coefficients. The second step involves transferring the derivative 
from the basis function u to the testing function v via Green’s identity as 
follows 


L 


rV 2 u 


- / Vu • VvdS + / 
Jn Jr 


du 

v—dl 
ri+r 2 un> 


(25) 


Substituting the above identity into (24), we arrive at the weak form of the 
Helmholtz equation 

/ (Vu • Vu — k 2 uv\ dS = f v^dl (26) 

JQ •'Fi+r2 UTi 

where Tx and T 2 describe the boundary of the solution domain as shown in 
Figure 1. The absorbing boundary condition for two dimensions is to be 
applied at the contour boundary T 2 ) and must be expressed in terms of the 
normal derivative of the wave function. For a Galerkin formulation of (26), 
we replace the testing function v with the basis function u. 


3.2 Three-dimensional case 

For simplicity, we consider the problem of 3-D scattering by a conducting 
body. To solve this problem using the finite element method, it is necessary 
to enclose the scatterer within a fictitious surface denoted by S. Since the 
ABCs are derived only for spherical surfaces, the fictitious outer surface in 
our case is usually a sphere. The scattered field, denoted by E, satisfies 
the Helmholtz vector equation interior to S and the absorbing boundary 
condition at S. The equivalent variational problem for this is given by 

SF( E) = 0 
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where F denotes the functional given by [6] 

F(E) = ^ [(V x E) • (V x E) — k 2 E ■ E] dV 

, + J s [<*E? + /?( V x E) r 2 - /3(V • E t ) 2 ] dS (27) 

in which V denotes the volume enclosed by S. To discretize this functional, 
the volume V is subdivided into a number of small tetrahedra, each occupying 
the volume V e (e = 1,2 , •••,M), where M denotes the total number of 
elements. Within each element, the electric field is expressed as 

m 

E e = Y, E j W j = {W e } T {E e } = {£ e } T {W e } (28) 

j=i 

where Wy are the edge-based vector basis functions given in [2], Ej denote the 
expansion coefficients of the basis, m represents the number of edges in the 
element and the superscript stands for the element number. On substituting 
(28) into (27), we obtain 

M Ms 

F = £{£'}V]{£'} + £{E'} T [B']{£ ! } (29) 

e=l s=l 

where M a denotes the number of triangular surface elements on S. Also, the 
elements of the matrices [ A e ] and [i? s ] are given by 

A h= L [(V X W?) . (V x W*) - klY/lWpV 

Bti= f [«WJ • W-, + /J(V x WJ) r • (V x W’)„ - «V • W?,)(V ■ WJ, )]<iS 

«/ Og 

where V e denotes the volume of the eth tetrahedron and S a denotes the 
surface area of the sth triangular surface element on S. On assembling all 
M elements, (29) can also be written as 

F = {E} t [I<){E} 

where [K] is a N x N sparse, symmetric matrix with N being the total 
number of edges in the geometry and {£}isaiVxl column vector denoting 
the edge fields. The system of equations is then obtained by applying the 
Rayleigh-Ritz procedure which amounts to differentiating F with respect to 
each edge field and then setting the differential to zero. This yields 

(/<•]{£} = (0) (30) 

Upon imposing the boundary condition at the surface of the scatterer which 
states that the total tangential electric field vanishes at the conducting sur- 
face, we obtain the final system whose solution yields the field components 
over the entire domain. In our implementation, only the non-zero elements 
of the \K\ matrix were stored and the resultant system was solved using the 
bi conjugate gradient algorithm. 
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4 Results 


Finite element codes have been extensively tested using two-dimensional ab- 
sorbing boundary conditions and have been known to perform quite well for 
circular mesh termination boundaries. In figure 2(a), we consider a cylindri- 
cal problem[7] to compare the BGT and Engquist-Majda ABCs. This prob- 
lem is interesting since strong resonant phenomena are observed for small 
bodies. It consists of a conducting cylinder with a radius of one free-space 
wavelength covered by a half free-space wavelength thick dielectric with a rel- 
ative permittivity of 2 and a relative permeability of 2. Figure 2(b) shows the 
comparison of the bistatic pattern (0, nc = 180°) of two solutions of the BGT 
ABCs, with the outer boundary at a distance of 1.5 A and .5 A away from the 
scatterer. As expected, the 1.5 A distant boundary produces better results 
than the .5A distant one. Figure 2(c) plots the performance of the Engquist- 
Majda ABCs for the same geometry and angle of incidence. The BGT ABC 
with its outer boundary a mere .5A away outperforms the Engquist-Majda 
ABC when its outer boundary is 2.5A distant from the scatterer. 

Figure 3(a) analyses two conducting wedges, 5A 0 long, with the area be- 
tween them filled with a dielectric having c r = 2 and /x r = 2. The incident 
angle of the incoming plane wave is 0°. Again, the results obtained using the 
BGT absorbing boundary condition in figure 3(b) matches the data obtained 
from a hybrid finite element / boundary element analysis. Figure 3(c) shows 
the corresponding results obtained with the Engquist-Majda ABCs. 

In figure 4, we present the backscatter data from a 0.3A x 0.3A plate. 
The Sommerfeld radiation condition was employed at the mesh termination 
boundary placed 0.3A away from the edge. The comparison with a standard 
integral equation code is seen to be excellent. The backscatter pattern of 
the same plate is shown in figure 5 with a conformal outer boundary. The 
zero thickness plate lying at z — 0 on the x — y plane was enclosed within 
' a boundary resembling a cut-off sphere. This was achieved by making a 
plane cut at z — ±.3A of a .45A radius sphere and setting r in (23) to oo 
over the planar surfaces. Storage was reduced by about 10% on using the 
conformal outer boundary but the reduction should be more striking for 
larger geometries. 

Figure 6 compares the measured bistatic cross-section (0,- = 180°,</>,- = 
90°) of a metallic cube having an edge length of 0.755A with the correspond- 
ing pattern computed by the three-dimensional FE-ABC code. The second- 
order vector ABC was employed at the mesh truncation boundary which was 
placed only 0.1 A from the edge of the cube. About 33000 unknowns were 
used for the discretization of the computational domain and the [A] matrix 
contained a total of 264000 distinct non-zerc ,-ntries. The storage require- 
ment of this matrix is consequently much smaller than the 2000 unknown 
system generated by the boundary integral approach which has 4 million 
non-zero entries. 

Figure 7 presents backscatter data for a cylinder of radius 0.3 A and height 
0.6A. The data from the three-dimensional finite element code again compare 
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well with that obtained from a moment method-body of revolution code. The 
mesh was terminated at a distance of 0.3A from the edge of the scatterer and 
the system consisted of about 33000 unknowns and converged to the solution 
in about 350 iterations when using the Sommerfeld radiation condition. Each 
iteration took approximately 0.3 seconds on a Cray YMP and on the average 
it was found that for N > 25000, the number of iterations required was of 
the order of N/ 100. 

Figure 8 shows the backscatter pattern on enclosing the pec cylinder 
in figure 7 with a conformal cylindrical boundary capped at the ends by a 
spherical surface. The storage was reduced by about 20% but the computed 
values are off by about ldB for most incidence angles. The poor results 
are probably due to the fact that the ABCs were derived only for spherical 
surfaces and therefore fail on cylindrical boundaries. 

The above results point out that the mesh truncation boundary could 
be placed a few tenths of a wavelength away from the scatterer to obtain 
reliable results. This is probably because in three dimensions, the scattered 
fields decrease as 1/r whereas in two dimensions, the rate of decrease is 
proportional to l/y/r. 
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Figure 1: Geometry of p.e.c scatterer enclosed in an outer circular boundary 




(a) 



(b) 



(c) 


Figure 2: (a)Cylinder geometry (b)Bistatic pattern using BGT absorbing 
boundary conditions (c)Bistatic pattern using Engquist-Majda absorbing 
boundary conditions 
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Figure 3: (a)Wedge geometry (b)Bistatic pattern using BGT absorbing 
boundary conditions (c)Bistatic pattern using Engquist-Majda absorbing 
boundary conditions 
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Figure 5: Backscatter pattern of the plate in figure 4 enclosed by a conformal 
outer boundary and using the 2nd order absorbing boundary condition. 
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Figure 6: Bistatic pattern (0; = 180°) of a metallic cube having an edge 
length of 0.755A using the 2nd order vector absorbing boundary condition. 
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Figure 7: Backscatter pattern of a metallic cylinder of radius 0.3A and height 
0.6A using the Sommerfeld radiation condition. The axis of the cylinder 
coincides with the z-ajcis of* the coordinate system. 
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Figure 8: Backscatter pattern of the cylinder in figure 7 enclosed by a con- 
formal outer boundary and using the second order absorbing boundary con- 
dition. 



